Effects of Abortion Bans

options(dplyr.summarise.inform = FALSE)

library(tidyverse)
library(tidybayes)
library(posterior)
library(jsonlite)
library(kableExtra)
library(gt)

suffix <- "main_analysis"
fig_suffix <- ""

df <- read_csv(sprintf("results/df_%s.csv", suffix))
aggregations <- read_csv("data/dobbs_biannual_aggragtions.csv")

df$time <- paste(df$year, df$bacode * 6 - 5, "01", sep="-") %>% ymd()
aggregations$time <- paste(aggregations$year, aggregations$bacode * 6 - 5, "01", sep="-") %>% ymd()
df <- df %>% filter(time >= "2012-01-01")
aggregations <- aggregations %>% filter(time >= "2012-01-01")

df |>
  mutate(
          start_date = ym(paste(year, "-", bacode * 6 - 5)),
          end_date = start_date + months(6) - days(1)
        ) -> df


fill_in_missing_denoms <- function(dat) {
    pop_index_2022 <- which.max(dat$year == 2022)
    pop_index_2021 <- which.max(dat$year == 2021)
    dat %>% mutate_at(vars(contains("pop")), ~ ifelse(is.na(.), .[pop_index_2022]^2 / .[pop_index_2021], .))
}

## Hacky imputation
df <- df %>%
    group_by(state) %>%
    group_modify(~ fill_in_missing_denoms(.)) %>%
    ungroup()

df$dobbs_code <- df$dobbscodev2
df <- df %>% group_by(state) %>% fill(exposed_infdeaths, .direction="down") %>% ungroup()

df <- df %>% group_by(state) %>% fill(exposed_infdeaths, .direction="down") %>% ungroup()
df <- df %>% group_by(state) %>% mutate(ban = ifelse(any(exposed_infdeaths == 1, na.rm=TRUE), TRUE, FALSE))


df$births_con <- df$births_noncon <- df$births_total

df$births_neo <- df$births_nonneo <- df$births_total
df$deaths_nonneo <- df$deaths_total - df$deaths_neo

agg_births <- df %>% select(state, year, bacode, time, starts_with("births"), exposed_infdeaths) %>% 
    group_by(state) %>%
    mutate(ban = ifelse(any(exposed_infdeaths == 1), "Exposed (excl. Texas)", "Unexposed")) %>%
    mutate(ban = ifelse(is.na(ban), "Unexposed", ban)) %>% 
    mutate(ban = ifelse(state == "Texas", "Texas", ban)) %>% 
    ungroup() %>% 
    pivot_longer(starts_with("births"), names_pattern="births_(.*)", names_to="category", values_to="births") %>%
    group_by(ban, time, category) %>% 
    summarize(births=sum(births)) %>% ungroup() %>% 
    mutate(state = ban)

agg_deaths <- df %>% select(-starts_with("births")) %>% 
    filter(state == "Texas") %>% bind_rows(aggregations %>% mutate(state = expcat)) %>% 
    pivot_longer(cols=starts_with("deaths"), names_pattern = "deaths_(.*)", names_to="category", values_to="deaths") %>%
    select(time, category, deaths, state, exposed_infdeaths) %>%
    mutate(state = recode(state, "exp"="Exposed (excl. Texas)", "Texas"="Texas", "unexp"="Unexposed"))

agg_df <- left_join(agg_deaths, agg_births, by=c("state", "time", "category"))

agg_df_total <- agg_df %>% mutate(mortality_rate = deaths / births * 1000) %>%
  filter(category == "total") %>%
  mutate(ban = (time >= "2022-07-01") & (state == "Texas") | 
               (time >= "2023-01-01") & (state == "Exposed (excl. Texas)"))



agg_df_total <- agg_df_total %>% arrange(time)


## span = 0.5
smooth_texas <- predict(loess(mortality_rate ~ as.numeric(time), data = agg_df_total %>% filter(state == "Texas") %>% arrange(time), span=0.5))
smooth_unexposed <- predict(loess(mortality_rate ~ as.numeric(time), data = agg_df_total %>% filter(state == "Unexposed") %>% arrange(time), span=0.5))
smooth_exposed <- predict(loess(mortality_rate ~ as.numeric(time), data = agg_df_total %>% filter(state == "Exposed (excl. Texas)") %>% arrange(time), span=0.5))


agg_df_total$smooth <- numeric(nrow(agg_df_total))

agg_df_total$smooth[agg_df_total$state == "Texas"] <- smooth_texas
agg_df_total$smooth[agg_df_total$state == "Exposed (excl. Texas)"] <- smooth_exposed
agg_df_total$smooth[agg_df_total$state == "Unexposed"] <- smooth_unexposed


tmp_tx <- agg_df_total %>% filter(state == "Texas", time == "2022-07-01")
tmp_tx$ban = FALSE
tmp_exp <- agg_df_total %>% filter(state == "Exposed (excl. Texas)", time == "2023-01-01")
tmp_exp$ban = FALSE
agg_df_total <- bind_rows(agg_df_total, tmp_tx, tmp_exp)
agg_df_total <- agg_df_total %>% arrange(time)

agg_df_total %>%
  ggplot() + 
  geom_smooth(aes(x=time, y=mortality_rate, col=state), span=0.5, se=FALSE) + 
  geom_line(aes(x=time, y=mortality_rate, col=state), alpha=0.5) +
  theme_bw(base_size=16) + 
  theme(legend.position = c(0.99, 0.99), legend.justification = c(1, 1),
  #legend.background = element_blank(),  # Make legend background transparent
  legend.title = element_blank()  ) +
  scale_color_manual(values=c("red", "orange", "dark gray"),
                      labels=c("States with bans (excl. Texas)", "Texas", "States without bans")) + 
                
  ylab("Infant Mortality Rate (per 1000 live births)") + xlab("Year") +
  scale_x_date(
      date_breaks = "1 year",  # Set axis labels at yearly intervals
      date_labels = "%Y",
      limits=c(as.Date("2012-01-01"), as.Date("2023-08-01"))
  ) + guides(linetype="none") +
  geom_vline(xintercept=lubridate::date("2022-01-01"), color="orange", linetype="dashed") +
  geom_vline(xintercept=lubridate::date("2023-01-01"), color="red", linetype="dashed") +
  geom_text(x=lubridate::date("2022-02-01"), y=6.2, label="Texas exposure in effect", vjust=1, 
  angle=90, color="orange", size=2.5) + 
  geom_text(x=lubridate::date("2023-02-01"), y=6.4, label="Other exposures in effect", vjust=1, 
  angle=90, color="red", size = 2.5) +
  coord_cartesian(xlim=c(as.Date("2012-04-01"), as.Date("2023-04-01")))

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

ggsave(sprintf("figs/%smain_figures/absolute_mortality.png", fig_suffix), width=8, height=5)



relative_mortality_rate <- agg_df %>% 
    mutate(mortality_rate = deaths / births * 1000) %>%
    group_by(state, category) %>% 
    mutate(mean_mr = mean(mortality_rate[time < "2022-03-01"])) %>%
    mutate(relative_mr = mortality_rate / mean_mr)

relative_mortality_rate %>% filter(category == "total") %>% 
    ggplot() + 
    geom_smooth(aes(x=time, y=relative_mr, group=state, col=state), se=FALSE) + 
    geom_point(aes(x=time, y=relative_mr, col=state), height=0, alpha=0.5) +
    theme_bw(base_size=16) + 
    theme(legend.position = c(0.99, 0.99), legend.justification = c(1, 1),
    #legend.background = element_blank(),  # Make legend background transparent
    legend.title = element_blank()  ) +
    scale_color_manual(values=c("red", "orange", "dark gray"),
                       labels=c("Banned States (excl. Texas)", "Texas", "Non-banned States")) + 
    ylab("Relative Infant Mortality Rate") + xlab("Year") + 
    geom_vline(xintercept=lubridate::date("2022-01-01"), color="orange", linetype="dashed") +
    geom_vline(xintercept=lubridate::date("2023-01-01"), color="red", linetype="dashed") +
    scale_x_date(
        date_breaks = "1 year",  # Set axis labels at yearly intervals
        date_labels = "%Y"       # Format the labels to show the year
    ) 

Warning in geom_point(aes(x = time, y = relative_mr, col = state), height = 0,
: Ignoring unknown parameters: `height`

ggsave(sprintf("figs/%smain_figures/relative_mortality.png", fig_suffix), width=8, height=5)

relative_mortality_rate %>%
    ggplot() + 
    geom_smooth(aes(x=time, y=relative_mr, group=state, col=state), se=FALSE) + 
    geom_point(aes(x=time, y=relative_mr, col=state), height=0, alpha=0.5) +
    theme_bw(base_size=16) + 
    theme(legend.position = c(0.99, 0.99), legend.justification = c(1, 1),
    #legend.background = element_blank(),  # Make legend background transparent
    legend.title = element_blank()  ) +
    scale_color_manual(values=c("red", "orange", "dark gray")) + 
    ylab("Relative Mortality Rate") + xlab("Year") + 
    geom_vline(xintercept=lubridate::date("2022-01-01"), color="orange", linetype="dashed") +
    geom_vline(xintercept=lubridate::date("2023-01-01"), color="red", linetype="dashed") +
    facet_wrap(~category, scales="free_y")

Warning in geom_point(aes(x = time, y = relative_mr, col = state), height = 0,
: Ignoring unknown parameters: `height`

ggsave(sprintf("figs/%smain_figures/relative_facet.png", fig_suffix), width=14, height=10)

categories_list <- list(
    race = c("nhwhite", "nhblack", "hisp", "otherraceeth"),
    total = c("total"),
    neonatal = c("neo", "nonneo"),
    congenital = c("con", "noncon")
)

########################
## Missing Data Table
########################

df %>% select(time, state, starts_with("deaths"), ban) %>% 
  filter(ban == 0) %>%                   
  pivot_longer(cols=starts_with('deaths'), names_pattern="deaths_(.*)", 
                    names_to="category", 
                    values_to="deaths") %>%
  group_by(category) %>% 
  summarize(count_missing = sum(is.na(deaths)), fraction_missing = round(100*count_missing / n(), 1)) %>%
  mutate(state = "States without bans") -> frac_missing_unexposed

df %>% select(time, state, starts_with("deaths"), ban) %>% 
  filter(ban == 1 & state != "Texas") %>%                   
  pivot_longer(cols=starts_with('deaths'), names_pattern="deaths_(.*)", 
                    names_to="category", 
                    values_to="deaths") %>%
  group_by(category) %>% 
  summarize(count_missing = sum(is.na(deaths)), fraction_missing = round(100*count_missing / n(), 1)) %>%
  mutate(state = "States with bans (excl. Texas)") -> frac_missing_ban_no_tx

df %>% select(time, state, starts_with("deaths"), ban) %>% 
  filter(ban == 1) %>%                   
  pivot_longer(cols=starts_with('deaths'), names_pattern="deaths_(.*)", 
                    names_to="category", 
                    values_to="deaths") %>%
  group_by(category) %>% 
  summarize(count_missing = sum(is.na(deaths)), fraction_missing = round(100*count_missing / n(), 1)) %>% 
  mutate(state = "States with bans") -> frac_missing_ban_states

df %>% select(time, state, starts_with("deaths"), ban) %>% 
  filter(ban == 1) %>%                   
  pivot_longer(cols=starts_with('deaths'), names_pattern="deaths_(.*)", 
                    names_to="category", 
                    values_to="deaths") %>% 
   
  group_by(category, state) %>% 
  summarize(count_missing = sum(is.na(deaths)), fraction_missing = round(100*count_missing / n(), 1)) %>% 
  bind_rows(frac_missing_unexposed, frac_missing_ban_no_tx, frac_missing_ban_states) %>%
  select(category, state, fraction_missing) %>%
  pivot_wider(names_from = category, values_from = fraction_missing) |> 
  gt() |>
  tab_spanner(
    label = "Race",
    columns = c(nhwhite, nhblack, hisp, otherraceeth)) |>
  tab_spanner(
    label = "Timing",
    columns = c(neo, nonneo)) |>
    tab_spanner(
    label = "Type",
    columns = c(con, noncon)) |>
  cols_label(
   state = "",
   con = "Congenital",
   noncon = "Noncongenital",
   neo = "Neonatal",
   nonneo = "Nonneonatal",
   nhwhite = "Non-Hispanic White",
   nhblack = "Non-Hispanic Black",
   hisp = "Hispanic",
   otherraceeth = "Other",
   total = "Total"
   ) |>
  tab_row_group(
    label = md("**States with bans**"),
    rows = 1:14,
    id = "states_w_bans"
  ) |>
  tab_row_group(
    label = md("**Aggregated**"),
    rows = 15:17,
    id = "aggregates"
  ) |>
   tab_style(style = cell_borders(sides = c("bottom"),  weight = px(1.5)),
            locations = cells_body(rows = c(14))) |>
   fmt_number(columns = total, decimals = 1) |>
   gtsave(sprintf("figs/%ssupplement_figures/missing_data.png", fig_suffix), zoom=4)


## Median infant death counts
agg_deaths %>% ungroup() %>%
  filter(state == "Unexposed") %>%             
  group_by(category) %>% 
  summarize(count_missing = round(mean(deaths, na.rm=TRUE), 0)) %>%
  mutate(state = "States without bans") -> median_deaths_unexposed

agg_deaths %>% ungroup() %>%
  filter(state == "Exposed (excl. Texas)") %>%             
  group_by(category) %>% 
  summarize(count_missing = round(mean(deaths, na.rm=TRUE), 0)) %>%
  mutate(state = "States with bans (excl. Texas)") -> median_deaths_ban_no_tx

agg_deaths %>% ungroup() %>%
  filter(state == "Exposed (excl. Texas)" | state == "Texas") %>%
  group_by(category, time) %>% 
  summarize(count_missing = sum(deaths, na.rm=TRUE)) %>%
  group_by(category) %>% 
  summarize(count_missing = round(mean(count_missing, na.rm=TRUE), 0)) %>%
  mutate(state = "States with bans") -> median_deaths_ban


df %>% select(time, state, starts_with("deaths"), ban) %>% 
  filter(ban == 1) %>%                   
  pivot_longer(cols=starts_with('deaths'), names_pattern="deaths_(.*)", 
                    names_to="category", 
                    values_to="deaths") %>% 
  group_by(category, state) %>% 
  mutate(deaths = replace_na(deaths, 5)) %>%
  summarize(count_missing = round(median(deaths, na.rm=TRUE), 0)) %>% 
  bind_rows(median_deaths_unexposed, median_deaths_ban_no_tx, median_deaths_ban) %>%
  mutate(count_missing = ifelse(count_missing < 10, "<10", as.character(count_missing))) %>%
  select(category, state, count_missing) %>%
  pivot_wider(names_from = category, values_from = count_missing) |> 
  gt() |>
  tab_spanner(
    label = "Race",
    columns = c(nhwhite, nhblack, hisp, otherraceeth)) |>
  tab_spanner(
    label = "Timing",
    columns = c(neo, nonneo)) |>
    tab_spanner(
    label = "Type",
    columns = c(con, noncon)) |>
  cols_label(
   state = "",
   con = "Congenital",
   noncon = "Noncongenital",
   neo = "Neonatal",
   nonneo = "Nonneonatal",
   nhwhite = "Non-Hispanic White",
   nhblack = "Non-Hispanic Black",
   hisp = "Hispanic",
   otherraceeth = "Other",
   total = "Total"
   ) |>
  tab_row_group(
    label = md("**States with bans**"),
    rows = 1:14,
    id = "states_w_bans"
  ) |>
  tab_row_group(
    label = md("**Aggregated**"),
    rows = 15:17,
    id = "aggregates"
  ) |>
   tab_style(style = cell_borders(sides = c("bottom"),  weight = px(1.5)),
            locations = cells_body(rows = c(14))) |>
  cols_align(
    align = "center",
    columns = -1
  ) |>
   gtsave(sprintf("figs/%ssupplement_figures/median_death_counts.png", fig_suffix), zoom=4)

source("plot_utilities.R")

## Load all datasets
file_dir <- "results"
types <- c("race", "total", "congenital", "neonatal")
ranks <- c(2, 2, 3, 3)

all_samples <- tibble()
for(i in 1:length(types)) {
    type <- types[i]
    model_rank <- ranks[i]
    print(sprintf("%s %i", type, model_rank))
    
    categories <- categories_list[[type]]
    type_samples <- read_csv(sprintf("%s/Poisson_mortality_%s_%i_%s.csv", file_dir, type, model_rank, suffix))

    ## Thinning for speeding things up, only for debugging
    ## type_samples <- type_samples %>% filter(.draw %% 10 == 0)

    categories <- categories_list[[type]]

    agg_category_name = "total"

    merged_df <- merge_draws_and_data(df, 
                                      type_samples, categories=categories,
                                      agg_category_name=agg_category_name, 
                                      outcome_string="deaths", 
                                      denom_string="births")

    df_ban_no_tx_treat_time <- merged_df %>% 
    filter(!(state %in% c("Texas", "Ban States")),  exposure_code==1) %>% pull(time) %>% min()
    df_ban_no_tx <- merged_df %>% filter(state == "Ban States" ) %>% 
    mutate(exposure_code = ifelse(time >= df_ban_no_tx_treat_time, 1, 0))
    df_tx <- merged_df %>% filter(state == "Texas")
    df_ban_no_tx <- df_ban_no_tx %>% 
    mutate(state  ="Ban States (excl. Texas)",
            ban = TRUE,
            exposure_code = ifelse(time >= df_ban_no_tx_treat_time, 1, 0),
            ypred = df_ban_no_tx$ypred - df_tx$ypred,
            mu = log(exp(df_ban_no_tx$mu) - exp(df_tx$mu)),
            mu_treated = log(exp(df_ban_no_tx$mu_treated) - exp(df_tx$mu_treated)),
            deaths = df_ban_no_tx$deaths - df_tx$deaths,
            births = df_ban_no_tx$births - df_tx$births
            )

    agg_long_exposed <- aggregations %>% pivot_longer(cols=starts_with("deaths"), names_to=c("category"), names_pattern="deaths_(.*)$", values_to="deaths") %>% filter(expcat == "exp")
    agg_long_unexposed <- aggregations %>% pivot_longer(cols=starts_with("deaths"), names_to=c("category"), names_pattern="deaths_(.*)$", values_to="deaths") %>% filter(expcat == "unexp")

    df_ban_no_tx <- left_join(df_ban_no_tx, agg_long_exposed, by=c("time"="time", "category"="category")) %>% mutate(deaths = deaths.y) 

    df_all_ban <- left_join(merged_df %>% filter(state=="Ban States") %>% select(-deaths), 
                            agg_long_exposed, by=c("time"="time", "category"="category")) %>% 
                  mutate(deaths = deaths + df_tx$deaths) %>% ## Texas Deaths + All Deaths excluding Texas
                  mutate(state = "Ban States",
                         ban = TRUE,
                         ypred = df_ban_no_tx$ypred + df_tx$ypred,
                         mu = log(exp(df_ban_no_tx$mu) + exp(df_tx$mu)),
                         mu_treated = log(exp(df_ban_no_tx$mu_treated) + exp(df_tx$mu_treated)))

    df_banned_states <- merged_df %>% filter(time < df_ban_no_tx_treat_time)

    df_unexposed <- merged_df %>% filter(state == "Unexposed States" ) %>% 
    mutate(exposure_code = 0)
    df_unexposed <- left_join(df_unexposed, agg_long_unexposed, by=c("time"="time", "category"="category")) %>% mutate(deaths = deaths.y)

    merged_df <- merged_df %>% filter(state != "Ban States") %>% bind_rows(df_all_ban %>% select(colnames(merged_df)))
    merged_df <- merged_df %>% bind_rows(df_ban_no_tx)
    merged_df <- merged_df %>% filter(state != "Unexposed States") %>% bind_rows(df_unexposed)

    merged_df$type <- type
    merged_df$rank <- model_rank
    all_samples <- bind_rows(all_samples, merged_df)
}

[1] "race 2"

New names:
Rows: 1000 Columns: 15347
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(15347): ...1, .chain, .draw, mu[1,1,1], mu[2,1,1], mu[3,1,1], mu[4,1,1]...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

[1] "total 2"

New names:
Rows: 1000 Columns: 3878
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(3878): ...1, .chain, .draw, mu[1,1,1], mu[1,2,1], mu[1,3,1], mu[1,4,1],...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

[1] "congenital 3"

New names:
Rows: 1000 Columns: 7851
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(7851): ...1, .chain, .draw, mu[1,1,1], mu[2,1,1], mu[1,2,1], mu[2,2,1],...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

[1] "neonatal 3"

New names:
Rows: 1000 Columns: 7851
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(7851): ...1, .chain, .draw, mu[1,1,1], mu[2,1,1], mu[1,2,1], mu[2,2,1],...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

all_samples <- all_samples %>% mutate(category = fct_recode(category,
                                "Hispanic" = "hisp",
                                "Non-Hispanic Black" = "nhblack",
                                "Non-Hispanic White" = "nhwhite",
                                "Other" = "otherraceeth",
                                "Total" = "total",
                                "Congenital" = "con",
                                "Non-congenital" =  "noncon",
                                "Neonatal" = "neo",
                                "Non-neonatal" = "nonneo",
                                )) %>% 
                                mutate(category = fct_relevel(category,
                                "Congenital", "Non-congenital", "Neonatal", "Non-neonatal",
                                "Hispanic", "Non-Hispanic Black", "Non-Hispanic White", "Other",
                                "Total"
                                ))

quantiles_df <- all_samples %>% group_by(category, type, state, time) %>%
  summarize(ypred_mean=mean(ypred), 
            ypred_lower=quantile(ypred, 0.025), ypred_upper=quantile(ypred, 0.975), 
            births=mean(births), 
            deaths=mean(deaths), 
            exposure_code = first(exposure_code),
            ban = first(ban)) %>% ungroup()

Fit and Gap Plots

make_state_fit_plot(quantiles_df %>% filter(type=="total"), state_name="Ban States", category="Total", target="deaths") + theme_bw(base_size=16) +  
scale_x_date(
    date_breaks = "1 year",  # Set axis labels at yearly intervals
        date_labels = "%Y"       # Format the labels to show the year
) +   
theme(
    plot.margin = margin(t = 10, r = 20, b = 10, l = 10)  # Extend the right margin
) + ggtitle("Model Fit - All states with bans")

ggsave(sprintf("figs/%smain_figures/ban_states_fit_plot.png", fig_suffix), width=8, height=5)

make_gap_plot(quantiles_df %>% filter(type=="total"), state_name="Ban States", category="Total", target="deaths") + theme_bw(base_size=16) + 
scale_x_date(
    date_breaks = "1 year",  # Set axis labels at yearly intervals
    date_labels = "%Y"       # Format the labels to show the year
) +   
ggtitle("Gap Plot - All states with bans") +
theme(
    plot.margin = margin(t = 10, r = 20, b = 10, l = 10)  # Extend the right margin
)

ggsave(sprintf("figs/%smain_figures/ban_states_gap_plot.png", fig_suffix),width=8, height=5)

Interval Plots

make_interval_plot(all_samples %>% filter(type == "race", category !="Total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu") + labs(color="Race and ethnicity") +
                         facet_grid(factor(!(state %in% c("States w/ bans", "States w/ bans (excl. Texas)"))) ~ . , scales="free_y", space='free') + ggtitle("Race and ethnicity")

ggsave(sprintf("figs/%smain_figures/race_interval_plot.png", fig_suffix), width=10, height=10)

race_interval_agg <- make_interval_plot(
  all_samples %>% filter(type == "race", category !="Total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu",
                         states = c("Ban States", "Ban States (excl. Texas)", "Texas")) + 
                         labs(color="Race and ethnicity") + ggtitle("Race and ethnicity") + ylim(c(-5, 25))


make_interval_plot(all_samples %>% filter(type == "congenital", category !="Total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu") + labs(color="Congenital") + 
                         facet_grid(factor(!(state %in% c("States w/ bans", "States w/ bans (excl. Texas)"))) ~ . , scales="free_y", space='free')

ggsave(sprintf("figs/%smain_figures/congenital_interval_plot.png", fig_suffix), width=10, height=10)

congenital_interval_agg <- make_interval_plot(all_samples %>% filter(type == "congenital", category !="Total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu",
                         states = c("Ban States", "Ban States (excl. Texas)", "Texas")) + labs(color="Congenital") + 
                         ggtitle("Type of Death") + ylim(limits=c(-5, 25))


make_interval_plot(all_samples %>% filter(type == "neonatal", category !="Total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu") + labs(color="Neonatal") + 
                         facet_grid(factor(!(state %in% c("States w/ bans", "States w/ bans (excl. Texas)"))) ~ . , scales="free_y", space='free')

ggsave(sprintf("figs/%smain_figures/neonatal_interval_plot.png", fig_suffix), width=10, height=10)

neonatal_interval_agg <- 
make_interval_plot(all_samples %>% filter(type == "neonatal", category !="Total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu", 
                         states = c("Ban States", "Ban States (excl. Texas)", "Texas")) + labs(color="Neonatal") + 
                         ggtitle("Timing of Death") + ylim(c(-5, 25))



make_interval_plot(all_samples %>% filter(type == "total"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu") + labs(color="Total") + 
                         facet_grid(factor(!(state %in% c("States w/ bans", "States w/ bans (excl. Texas)"))) ~ . , scales="free_y", space='free') + ggtitle("All Infant Deaths")

ggsave(sprintf("figs/%smain_figures/total_interval_plot.png", fig_suffix), width=10, height=10)


make_interval_plot(all_samples %>% filter(type == "race", category =="Non-Hispanic Black"), 
                         group_var = c("state", "category"),
                         target="deaths", denom="births",
                         rate_normalizer=1000,
                         estimand = "ratio", 
                         method="mu") + labs(color="Race and ethnicity") +
                         facet_grid(factor(!(state %in% c("States w/ bans", "States w/ bans (excl. Texas)"))) ~ . , scales="free_y", space='free') + ggtitle("Non-Hispanic Black Infant Deaths")

ggsave(sprintf("figs/%smain_figures/black_interval_plot.png", fig_suffix), width=10, height=10)

## All categories for banned states
## Interval Plot For Table 1
make_interval_plot(all_samples %>% 
                   filter((type == "total" & category == "Total") | category != "Total") %>% 
                   mutate(category = fct_rev(category), type=fct_rev(type)), 
                   states = c("Ban States"),
                   group_var = c("type", "category", "state"),
                   target="deaths", denom="births",
                   rate_normalizer=1000,
                   estimand = "ratio", 
                   method="mu", 
                   color_group = "state",
                   x_var = "category") + 
                   labs(color="Total") + 
                   facet_grid(type ~ ., scales="free_y", space='free') +
                   theme(panel.spacing = unit(2, "lines"))

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `state = fct_relevel(state, "Ban States (excl. Texas)", "Ban
  States")`.
Caused by warning:
! 1 unknown level in `f`: Ban States (excl. Texas)

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `state = fct_recode(state, `States w/ bans` = "Ban States",
  `States w/ bans (excl. Texas)` = "Ban States (excl. Texas)")`.
Caused by warning:
! Unknown levels in `f`: Ban States (excl. Texas)

ggsave(sprintf("figs/%smain_figures/ban_states_interval_plot.png", fig_suffix), width=10, height=6)

race_interval_agg / congenital_interval_agg / neonatal_interval_agg

Warning: Removed 183 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Warning: Removed 112 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Warning: Removed 6 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

ggsave(sprintf("figs/%smain_figures/agg_interval_plot.png", fig_suffix), width=10, height=10)

Warning: Removed 183 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Warning: Removed 112 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Warning: Removed 6 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Table

make_mortality_table <- function(merged_df, 
                       target_state = "Texas", target="deaths", denom="births",
                       rate_normalizer=1000, 
                       tab_caption = "Table 1. Expected difference† in infant deaths (count and rate) in all states that banned abortion in months affected by bans, overall, by race, timing of death, and cause of death.") {

  if(target_state == "Ban States") {
     
     ###################
     ## Before 2023 Ban States is just Texas counts
     ###################
     merged_df_tx <- merged_df %>% filter(state == "Texas")
     merged_df <- merged_df %>% filter(state == "Ban States")     
     merged_df[merged_df$time < "2023-01-01", c("deaths", "births", "ypred", "mu", "mu_treated")] <- merged_df_tx[merged_df_tx$time < "2023-01-01", c("deaths", "births", "ypred", "mu", "mu_treated")]
     merged_df <- merged_df %>% filter(exposure_code == 1) %>% mutate(denom = births)
    #  merged_df <- merged_df %>%
    #    filter(exposure_code == 1) %>%
    #    ## Aggregate over all banned states
    #    group_by(type, category, .draw, time) %>% 
    #    summarise({{target}} := sum(.data[[target]]), 
    #              denom = sum(.data[[denom]]), 
    #              ypred=sum(ypred), 
    #              mu = log(sum(exp(mu))),
    #              mu_treated = log(sum(exp(mu_treated))),
    #              years=mean(interval(start_date, end_date) / years(1)))
  } else {
    merged_df <- merged_df %>%
      filter(state == target_state, exposure_code == 1) %>%
      mutate(denom=.data[[denom]])
  }  
#   } else if(target_state == "Ban States (excl.Texas)") {
#     merged_df <- merged_df %>% filter(!state %in% c("Ban States", "Ban States (excl. Texas)"))
#     merged_df <- merged_df %>%
#       filter(state != "Texas") %>%
#       filter(exposure_code == 1) %>%
#       ## Aggregate over all banned states
#       group_by(type, category, .draw, time) %>% 
#       summarise({{target}} := sum(.data[[target]]), 
#                 denom = sum(.data[[denom]]), 
#                 ypred=sum(ypred), 
#                 mu = log(sum(exp(mu))),
#                 mu_treated = log(sum(exp(mu_treated))),
#                 years=mean(interval(start_date, end_date) / years(1)))


  
  table_df <- merged_df %>%
    ungroup() %>%
    ## Aggregate over time
    group_by(type, category, .draw) %>%
    summarize(
      ypred = sum(ypred),
      outcome = sum(ifelse(is.na(.data[[target]]), round(exp(mu)), .data[[target]])), 
      treated = sum(exp(mu_treated)), untreated = sum(exp(mu)),
      denom = sum(denom, na.rm = TRUE),
      treated_rate = treated / denom * rate_normalizer,
      untreated_rate = untreated / denom * rate_normalizer,
      outcome_rate = round(outcome / denom * rate_normalizer, 2),
      outcome_diff = round(treated - untreated)
    ) %>%
    ungroup() %>%
    ## Compute quantiles of effects
    group_by(type, category) %>%
    summarize(
      ypred_mean = mean(ypred),
      outcome = round(mean(outcome), digits=ifelse(mean(outcome < 10), 2, 0)),
      outcome_diff_mean = round(mean(outcome_diff), digits=ifelse(mean(outcome < 1), 2, 0)), 
      outcome_diff_lower = round(quantile(outcome_diff, 0.025)), 
      outcome_diff_upper = round(quantile(outcome_diff, 0.975)),
      outcome_rate = mean(outcome_rate),
      ypred_lower = quantile(ypred, 0.025), ypred_upper = quantile(ypred, 0.975),
      treated_mean = mean(treated), treated_lower = quantile(treated, 0.025), treated_upper = quantile(treated, 0.975),
      untreated_mean = mean(untreated), untreated_lower = quantile(untreated, 0.025), untreated_upper = quantile(untreated, 0.975),      
      treated_rate_mean = mean(treated_rate), treated_rate_lower = quantile(treated_rate, 0.025), treated_rate_upper = quantile(treated_rate, 0.975),
      untreated_rate_mean = mean(untreated_rate), untreated_rate_lower = quantile(untreated_rate, 0.025), untreated_rate_upper = quantile(untreated_rate, 0.975), 
      causal_effect_diff_mean = mean(treated_rate - untreated_rate), causal_effect_diff_lower = quantile(treated_rate - untreated_rate, 0.025), causal_effect_diff_upper = quantile(treated_rate - untreated_rate, 0.975),
      causal_effect_ratio_mean = mean(treated_rate / untreated_rate), causal_effect_ratio_lower = quantile(treated_rate / untreated_rate, 0.025), causal_effect_ratio_upper = quantile(treated_rate / untreated_rate, 0.975),
      denom = mean(denom),
      pval = 2*min(mean(untreated_rate > treated_rate), mean(untreated < treated))
     )
    
  table_df <- table_df %>%
  mutate(
    # ypred_mean_rate = ypred_mean / years / (denom / rate_normalizer),
    outcome_rate = round(outcome_rate, 2),
    rate_diff = round(causal_effect_diff_mean, 2),
    rate_diff_lower = round(causal_effect_diff_lower, 2),
    rate_diff_upper = round(causal_effect_diff_upper, 2),
    mult_change = causal_effect_ratio_mean,
    mult_change_lower = causal_effect_ratio_lower,
    mult_change_upper = causal_effect_ratio_upper)

  table_df <- table_df %>%
    mutate(untreated_mean = round(untreated_mean, 0)) %>%
    mutate(untreated_rate_mean = round(untreated_rate_mean, 2)) %>% 
    mutate(death_counts_str = paste0(outcome_diff_mean, " (", outcome_diff_lower, ", ", outcome_diff_upper, ")")) %>%
    mutate(death_rate_abs_str = paste0(rate_diff, " (", rate_diff_lower, ", ", rate_diff_upper, ")")) %>%
    mutate(death_rate_pct_str = paste0(round(100*(mult_change-1), 2), " (", round(100*(mult_change_lower-1), 2), ", ", round(100*(mult_change_upper-1), 2), ")")) %>%
    ungroup() %>%
    filter((type != "total" & category !="Total")|type == "total")
  

  pvals <- pval_rows <- table_df %>% pull(pval)
  pval_rows <- which(pvals < 0.05)
  table_df <- table_df %>% mutate(category = paste0(category, ifelse(rate_diff_lower > 0, "*", "")))

  table_df %>%
    select(type, category, denom, outcome, outcome_diff_mean, outcome_rate, rate_diff, death_counts_str, death_rate_abs_str, death_rate_pct_str) %>%
    mutate(expected_outcome = outcome - outcome_diff_mean, expected_rate = outcome_rate - rate_diff) %>%
    mutate(outcome = as.character(outcome)) %>%
    select(-c("outcome_diff_mean", "rate_diff")) %>% 
    mutate(outcome = ifelse(as.numeric(outcome) < 10, NA, outcome)) %>%
    mutate(outcome_rate = ifelse(is.na(outcome), NA, outcome_rate)) %>%
    gt(rowname_col = "category") |>
  tab_header(
    title = tab_caption
  ) |> 
  ## ROW OPERATIONS
  tab_row_group(
    label = "Cause of death",
    rows = type == "congenital"
  ) |>
  tab_row_group(
    label = "Race and ethnicity",
    rows = type == "race"
  ) |>
  tab_row_group(
    label = "Timing",
    rows = type == "neonatal"
  ) |>
  row_group_order(groups = c(NA, "Race and ethnicity", "Timing",
                             "Cause of death")) |>
  ### COLUMN OPERATIONS
  tab_spanner(
    label = "Infant mortality rate (per 1,000 live births)",
    columns = c(outcome_rate, expected_rate, death_rate_abs_str, death_rate_pct_str)) |>
  tab_spanner(
    label = "Infant death count",
    columns = c(outcome, expected_outcome, death_counts_str)) |>
  cols_label(
    denom = "Births",
    outcome_rate = "Observed",
    expected_rate = "Expected",
    death_rate_abs_str = html("Expected difference<br>(95% CI)"),
    death_rate_pct_str = html("Expected percent change<br>(95% CI)"),
    outcome = "Observed",
    expected_outcome = "Expected",
    death_counts_str = html("Expected difference<br>(95% CI)"),
    category = ""
  ) |>
  tab_stub_indent(
    rows = category != "Total",
    indent = 5
  ) -> table_df
  
  ## Styling
  table_df |>
  tab_options(table.align = "left", heading.align = "left") |>
  cols_align(align = "left") |>
  cols_hide(c(type, category)) |>
  tab_options(table.font.size=8) |>
  opt_vertical_padding(scale = 0.5) |>
  cols_width(category ~ px(125),
             death_rate_abs_str ~ px(100),
             death_rate_pct_str ~ px(100),
             outcome_rate ~ px(50),
             death_counts_str ~ px(100),
             outcome ~ px(50),
             denom ~ px(50),
             expected_outcome ~ px(50),
             expected_rate ~ px(50)) |> 
  sub_missing() -> table_df_final

  table_df_final
}

ftab <- make_mortality_table(all_samples, target_state="Ban States") 
        # tab_footnote(html("Exposed months include January 2022 through December 2023 for Texas and January 2023 through December 2023 for other 13 states that banned abortion at 6 weeks or completely.<br>
        # † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.
        # The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br>
        # * Categories for which the 95% credible intervals of the expected difference excludes zero."))
ftab |> gtsave(sprintf("figs/%smain_figures/ban_states_table.png", fig_suffix), zoom=4)

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f52b9aa04f.html screenshot completed

ftab <- make_mortality_table(all_samples, target_state="Ban States (excl. Texas)", 
                             tab_caption = "Expected difference† in infant deaths (count and rate) in all states  excluding Texas that banned abortion in months affected by bans, overall, by race, timing of death, and cause of death.") 
        # tab_footnote(html("Exposed months include January 2023 through December 2023 for the 13 states besides Texas that banned abortion at 6 weeks or completely.<br>
        # † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.
        # The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br>
        # * Categories for which the 95% credible intervals of the expected difference excludes zero."))
ftab |> gtsave(sprintf("figs/%ssupplement_figures/ban_states_no_texas_table.png", fig_suffix), zoom=4)

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f53492b8e6.html screenshot completed

# ftab <- make_mortality_table(all_samples, target_state="Texas")
# ftab |> gtsave("~/Dropbox/abortion_results_and_data/mortality_figures/texas_table.png", zoom=4)

state_and_treat_times <- all_samples %>% filter(exposure_code==1, state != "Ban States") %>% 
                        group_by(state) %>% summarize(treatment_start = min(time))

treated_states <- state_and_treat_times %>% pull(state)
treated_times <- state_and_treat_times %>% pull(treatment_start)


for(i in 1:length(treated_states)) {
  target_state <- treated_states[i]

  if(target_state == "Ban States (excl. Texas)"){
    caption_string <- "all states with bans excluding Texas"
  }
  else{
    caption_string <- target_state
  }

  tab_caption = sprintf("Supplementary Table %i. Expected difference† in infant deaths (count and rate) in %s in months affected by bans, overall, by race, timing of death, and cause of death.", i, caption_string)
  
  make_mortality_table(all_samples, target_state=target_state, tab_caption=tab_caption) |>
  gtsave(filename=sprintf("figs/%ssupplement_figures/tables/%s_table.png", fig_suffix, target_state), zoom=4)


  # tab_footnote(html("† Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.
  #       The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br>
  #       * Categories for which the 95% credible intervals of the expected difference excludes zero. <br>
  #       In the publicly available CDC Wonder data, cell counts between 1 and 9 are suppressed. Dash denotes these missing counts or rates; expected differences are still inferred via the Bayesian hierarchical model.")) |>

  make_state_fit_plot(quantiles_df %>% filter(type=="total"), state_name=target_state, category="Total", target="deaths") + 
  theme_bw(base_size=16) + ggtitle(paste("Model Fit -", target_state, sep=" ")) + 
  scale_x_date(
    date_breaks = "1 year",  # Set axis labels at yearly intervals
        date_labels = "%Y"       # Format the labels to show the year
    ) +   
    theme(
        plot.margin = margin(t = 10, r = 20, b = 10, l = 10)  # Extend the right margin
    )
  ggsave(filename=sprintf("figs/%ssupplement_figures/fit_and_gap_plots/%s_fit_plot.png", fig_suffix, target_state), width=8, height=5)

  make_gap_plot(quantiles_df %>% filter(type=="total"), 
                state_name=target_state, category="Total", 
                target="deaths") + 
  theme_bw(base_size=16) + 
  ggtitle(paste("Gap Plot -", target_state, sep=" ")) +
  scale_x_date(
    date_breaks = "1 year",  # Set axis labels at yearly intervals
        date_labels = "%Y"       # Format the labels to show the year
    ) +   
    theme(
        plot.margin = margin(t = 10, r = 20, b = 10, l = 10)  # Extend the right margin
    )
  ggsave(sprintf("figs/%ssupplement_figures/fit_and_gap_plots/%s_gap_plot.png", fig_suffix, target_state), width=8, height=5)

}

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f5ca91f01.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f533328125.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f538bcbe65.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f56f77f5fc.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f52cb6ac90.html screenshot completed

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f5cfb38df.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f54152bf21.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f51f823108.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f51c611461.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f52980fad5.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f552d3cb3c.html screenshot completed

Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f54ec3fa9f.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f5247b051f.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f588d48ef.html screenshot completed

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f573af474b.html screenshot completed

States Table

make_state_table <- function(merged_df, 
                       target_category = "Total", target="deaths", denom="births",
                       rate_normalizer=1000, 
                       tab_caption = "Table %i. Expected difference† in infant deaths (count and rate) in states that banned abortion in months affected by bans.",
                       footnote_text="Exposed months include January 2022 through December 2023 for Texas and January 2023 through December 2023 for other 13 states that banned abortion at 6 weeks or completely.<br>
        † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.
        The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality). States are ordered by expected percent change.<br>
        * Categories for which the 95% credible intervals of the expected difference excludes zero.") {


  merged_df_tx <- merged_df %>% filter(state == "Texas")

  ## 2022-2023 "Texas is the only banned state, so all banned states is just Texas "
  ban_states_indices <- (merged_df$time < "2023-01-01" & merged_df$state == "Ban States")
  merged_df[ban_states_indices, c("deaths", "births", "ypred", "mu", "mu_treated")] <- 
  merged_df_tx[merged_df_tx$time < "2023-01-01", c("deaths", "births", "ypred", "mu", "mu_treated")]
  
  # ban_states_no_tx_indices <- (merged_df$time < "2023-01-01" & merged_df$state == "Ban States (excl. Texas)")
  # merged_df[ban_states_no_tx_indices, c("deaths", "births", "ypred", "mu", "mu_treated")] <-
  # (merged_df_tx[merged_df_tx$time >= "2023-01-01", c("deaths", "births", "ypred", "mu", "mu_treated")] - 
  # merged_df_tx[merged_df_tx$time >= "2023-01-01", c("deaths", "births", "ypred", "mu", "mu_treated")]


  merged_df <- merged_df %>% filter(exposure_code == 1, category == target_category) %>% mutate(denom = births)

  
  table_df <- merged_df %>%
    ungroup() %>%
    ## Aggregate over time
    group_by(state, .draw) %>%
    summarize(
      ypred = sum(ypred),
      outcome = sum(ifelse(is.na(.data[[target]]), round(exp(mu)), .data[[target]])), 
      treated = sum(exp(mu_treated)), untreated = sum(exp(mu)),
      denom = sum(denom, na.rm = TRUE),
      treated_rate = treated / denom * rate_normalizer,
      untreated_rate = untreated / denom * rate_normalizer,
      outcome_rate = round(outcome / denom * rate_normalizer, 2),
      outcome_diff = round(treated - untreated)
    ) %>%
    ungroup() %>%
    ## Compute quantiles of effects
    group_by(state) %>%
    summarize(
      ypred_mean = mean(ypred),
      outcome = round(mean(outcome), digits=ifelse(mean(outcome < 1), 2, 0)),
      outcome_diff_mean = round(mean(outcome_diff), digits=ifelse(mean(outcome < 1), 2, 0)), 
      outcome_diff_lower = round(quantile(outcome_diff, 0.025)), 
      outcome_diff_upper = round(quantile(outcome_diff, 0.975)),
      outcome_rate = mean(outcome_rate),
      ypred_lower = quantile(ypred, 0.025), ypred_upper = quantile(ypred, 0.975),
      treated_mean = mean(treated), treated_lower = quantile(treated, 0.025), treated_upper = quantile(treated, 0.975),
      untreated_mean = mean(untreated), untreated_lower = quantile(untreated, 0.025), untreated_upper = quantile(untreated, 0.975),      
      treated_rate_mean = mean(treated_rate), treated_rate_lower = quantile(treated_rate, 0.025), treated_rate_upper = quantile(treated_rate, 0.975),
      untreated_rate_mean = mean(untreated_rate), untreated_rate_lower = quantile(untreated_rate, 0.025), untreated_rate_upper = quantile(untreated_rate, 0.975), 
      causal_effect_diff_mean = mean(treated_rate - untreated_rate), causal_effect_diff_lower = quantile(treated_rate - untreated_rate, 0.025), causal_effect_diff_upper = quantile(treated_rate - untreated_rate, 0.975),
      causal_effect_ratio_mean = mean(treated_rate / untreated_rate), causal_effect_ratio_lower = quantile(treated_rate / untreated_rate, 0.025), causal_effect_ratio_upper = quantile(treated_rate / untreated_rate, 0.975),
      denom = mean(denom),
      pval = 2*min(mean(untreated_rate > treated_rate), mean(untreated < treated))
     )
    
  table_df <- table_df %>%
  mutate(
    # ypred_mean_rate = ypred_mean / years / (denom / rate_normalizer),
    outcome_rate = round(outcome_rate, 2),
    rate_diff = round(causal_effect_diff_mean, 2),
    rate_diff_lower = round(causal_effect_diff_lower, 2),
    rate_diff_upper = round(causal_effect_diff_upper, 2),
    mult_change = causal_effect_ratio_mean,
    mult_change_lower = causal_effect_ratio_lower,
    mult_change_upper = causal_effect_ratio_upper,
    state = fct_reorder(as.factor(state), mult_change, .fun = median),
    state = fct_recode(state, `States with bans` = "Ban States", `States with bans (excl. Texas)` = "Ban States (excl. Texas)")) %>%
    arrange(desc(state))
    
  

  table_df <- table_df %>%
    mutate(untreated_mean = round(untreated_mean, 0)) %>%
    mutate(untreated_rate_mean = round(untreated_rate_mean, 2)) %>% 
    mutate(death_counts_str = paste0(outcome_diff_mean, " (", outcome_diff_lower, ", ", outcome_diff_upper, ")")) %>%
    mutate(death_rate_abs_str = paste0(rate_diff, " (", rate_diff_lower, ", ", rate_diff_upper, ")")) %>%
    mutate(death_rate_pct_str = paste0(round(100*(mult_change-1), 2), " (", round(100*(mult_change_lower-1), 2), ", ", round(100*(mult_change_upper-1), 2), ")")) %>%
    ungroup()
  
  pvals <- pval_rows <- table_df %>% pull(pval)
  pval_rows <- which(pvals < 0.05)
  table_df <- table_df %>% mutate(state = paste0(state, ifelse(rate_diff_lower >= 0, "*", "")))
  
  table_df %>%
    select(state, denom, outcome, outcome_diff_mean, outcome_rate, rate_diff, death_counts_str, death_rate_abs_str, death_rate_pct_str) %>%
    mutate(expected_outcome = outcome - outcome_diff_mean, expected_rate = outcome_rate - rate_diff) %>%
    mutate(outcome = as.character(outcome)) %>%
    select(-c("outcome_diff_mean", "rate_diff")) %>% 
    gt(rowname_col = "state") |>
  tab_header(
    title = tab_caption
  ) |> 
  ## ROW OPERATIONS
  tab_row_group(
    label = md("**Aggregated**"),
    rows = str_detect(state, "States with bans")
  ) |>
  tab_row_group(
    label = md("**States with bans**"),
    rows = !str_detect(state, "States with bans")
  ) |>
  row_group_order(groups = c(md("**Aggregated**"), md("**States with bans**"))) |>

  ### COLUMN OPERATIONS
  tab_spanner(
    label = "Infant mortality rate (per 1,000 live births)",
    columns = c(outcome_rate, expected_rate, death_rate_abs_str, death_rate_pct_str)) |>
  tab_spanner(
    label = "Infant death count",
    columns = c(outcome, expected_outcome, death_counts_str)) |>
  cols_label(
    denom = "Births",
    outcome_rate = "Observed",
    expected_rate = "Expected",
    death_rate_abs_str = html("Expected difference<br>(95% CI)"),
    death_rate_pct_str = html("Expected percent change<br>(95% CI)"),
    outcome = "Observed",
    expected_outcome = "Expected",
    death_counts_str = html("Expected difference<br>(95% CI)"),
  ) -> table_df
  
  ## Styling
  table_df |>
  tab_options(table.align = "left", heading.align = "left") |>
  cols_align(align = "left") |>
  cols_hide(state) |>
  tab_options(table.font.size=8) |>
  opt_vertical_padding(scale = 0.5) |>
  cols_width(state ~ px(125),
             death_rate_abs_str ~ px(100),
             death_rate_pct_str ~ px(100),
             outcome_rate ~ px(50),
             death_counts_str ~ px(100),
             outcome ~ px(50),
             denom ~ px(50),
             expected_outcome ~ px(50),
             expected_rate ~ px(50)) -> table_df_final

  table_df_final |> tab_footnote(html(footnote_text))
}

total_states_tab <- make_state_table(all_samples %>% filter(type == "total"),
                    tab_caption = "Table 2. Expected difference† in total infant deaths (count and rate) split by states that banned abortion in months affected by bans.",
                    footnote_text = "")

# "Exposed months include January 2022 through December 2023 for Texas and January 2023 through December 2023 for other 13 states that banned abortion at 6 weeks or completely.<br>
#         † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.
#         The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br>
#         * Categories for which the 95% credible intervals of the expected difference excludes zero."
                    
total_states_tab |> gtsave(sprintf("figs/%smain_figures/states_table_total_deaths.png", fig_suffix), zoom=4)

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f5418048db.html screenshot completed

black_states_tab <- make_state_table(all_samples %>% filter(type == "race"),
                                    target_category="Non-Hispanic Black",
                                    tab_caption = "Table 3. Expected difference† in non-hispanic black infant deaths (count and rate) split by states that banned abortion in months affected by bans.",
                                    footnote_text = "")
# "Exposed months include January 2022 through December 2023 for Texas and January 2023 through December 2023 for other 13 states that banned abortion at 6 weeks or completely.<br>
#                                     † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.
#                                     The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br>
#                                     * Categories for which the 95% credible intervals of the expected difference excludes zero."
black_states_tab |> gtsave(sprintf("figs/%smain_figures/states_table_black_deaths.png", fig_suffix), zoom=4)

file:////var/folders/kg/96f7zf8d0wxdsgdpg521tsr00000gr/T//RtmpWFVNLk/file130f54d2f4375.html screenshot completed

PPC Figures

ppc_states <- c("Texas", "States with bans (excl. Texas)")
types <- c("race", "total", "congenital", "neonatal")
for(t in types) {

  figure_height = 10 * (length(categories_list[[t]])+1)/5

  if(t == "race") {
    all_samples_ppc <- all_samples %>% filter(category != "Other")
  } else {
    all_samples_ppc <- all_samples
  }

  all_samples_ppc <- all_samples_ppc %>% 
    mutate(state = fct_recode(state, `States with bans` = "Ban States", 
                                      `States with bans (excl. Texas)` = "Ban States (excl. Texas)")) 

  rmse_res <- make_rmse_ppc_plot(all_samples_ppc %>% filter(state %in% ppc_states, type == t), outcome="deaths")
  ggsave(rmse_res$rmse_plt, filename=sprintf("figs/%ssupplement_figures/ppc/%s_rmse_plot.png", fig_suffix, t), width=10, height=figure_height)

  abs_res <- make_abs_res_ppc_plot(all_samples_ppc %>% filter(state %in% ppc_states, type == t), outcome="deaths")
  ggsave(abs_res$max_plt, filename=sprintf("figs/%ssupplement_figures/ppc/%s_abs_plot.png", fig_suffix, t), width=10, height=figure_height)

  acf_ppc4 <- make_acf_ppc_plot(all_samples_ppc %>% filter(state %in% ppc_states, type == t),
                  lag=4, outcome="deaths") 
  ggsave(acf_ppc4$acf_plt, filename=sprintf("figs/%ssupplement_figures/ppc/%s_acf4_plot.png", fig_suffix, t), width=10, height=figure_height)

  acf_ppc2 <- make_acf_ppc_plot(all_samples_ppc %>% filter(state %in% ppc_states, type == t),
                  lag=2, outcome="deaths") 
  ggsave(acf_ppc2$acf_plt, filename=sprintf("figs/%ssupplement_figures/ppc/%s_acf2_plot.png", fig_suffix, t), width=10, height=figure_height)

  acf_ppc1 <- make_acf_ppc_plot(all_samples_ppc %>% filter(state %in% ppc_states, type == t), lag=1, outcome="deaths") 
  ggsave(acf_ppc1$acf_plt, filename=sprintf("figs/%ssupplement_figures/ppc/%s_acf1_plot.png", fig_suffix, t), width=10, height=figure_height)

  uc_ppcs_obj <- make_unit_corr_ppc_plot(all_samples_ppc %>% filter(type == t), outcome="deaths")
  ggsave(uc_ppcs_obj$eval_plt, filename=sprintf("figs/%ssupplement_figures/ppc/%s_unit_corr_plot.png", fig_suffix, t), width=7, height=figure_height)

}

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: There were 7992 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Hispanic` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [735] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 7991 remaining warnings.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: There were 1998 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Total` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [882] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1997 remaining warnings.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: There were 5994 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Congenital` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [882] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 5993 remaining warnings.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: There were 5994 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Neonatal` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [882] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 5993 remaining warnings.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Event-Time Gap Plot

all_samples %>% filter(ban==1, !(state %in% c("Ban States", "Ban States (excl. Texas)")), type=="total") %>% 
  mutate(event_time = as.numeric(factor(time, levels=unique(time)))) %>%
  group_by(state) %>% 
  mutate(min_event_time = min(event_time[exposure_code == 1])) %>%
  mutate(event_time = event_time - min_event_time) %>% 
  ungroup() %>% 
  group_by(event_time, .draw) %>% 
  summarize(deaths = sum(deaths), births=sum(births), 
            ypred = sum(ypred), mu = log(sum(exp(mu))), 
            mu_treated = log(sum(exp(mu_treated)))) %>% 
  ungroup() %>%
  group_by(event_time) %>%
  summarize(ypred_mean=mean(ypred), 
            ypred_lower=quantile(ypred, 0.025), ypred_upper=quantile(ypred, 0.975), 
            births=mean(births),
            deaths = mean(deaths)) %>% 
  ungroup() %>%
  mutate(pre_ban = event_time < 0) -> event_time_df

  event_time_df %>% write_csv("results/event_time_results.csv")

  event_time_df %>%
  ggplot() +
    geom_ribbon(aes(
      x = event_time,
      ymax = log(deaths) - log(ypred_lower),
      ymin = log(deaths) - log(ypred_upper),
      group = pre_ban
    ), alpha = 0.25) +
    geom_line(aes(x = event_time, y = log(deaths) - log(ypred_mean), group=pre_ban), color = "red") +
    theme_bw() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed", alpha=0.75) +
    xlab("Event Time")

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).

Aggregation Computations for Paper

file.create("figs/main_figures/summaries_for_paper.txt", overwrite=TRUE)

[1] TRUE TRUE

## Summaries based on fertility paper
mortality_expected_percent_change_total <- c(
  Texas = 9.41,
  Kentucky = 8.57,
  Alabama = 6.93,
  Oklahoma = 5.01,
  Arkansas = 3.22,
  Georgia = 2.01,
  Tennessee = 1.87,
  South_Dakota = 1.34,
  Mississippi = 1.31,
  Louisiana = -0.01,
  Wisconsin = -0.11,
  Idaho = -0.94,
  Missouri = -2.26,
  West_Virginia = -3.09
)

fertility_expected_percent_change_total <- c(
  Texas = 2.9,
  Kentucky = 1.64,
  Alabama = 1.27,
  Oklahoma = 1.36,
  Arkansas = 1.11,
  Georgia = 0.89,
  Tennessee = 1.11,
  South_Dakota = .63,
  Mississippi = 1.6,
  Louisiana = 1.2,
  Wisconsin = .85,
  Idaho = 0.73,
  Missouri = -.56,
  West_Virginia = .78
)

cor_test <- cor.test(mortality_expected_percent_change_total, fertility_expected_percent_change_total)

write_lines("------ Correlation between inferred pct increase in fertility with inferred pct increase in mortality  ------ ", "figs/main_figures/summaries_for_paper.txt", append=TRUE)
write_lines(sprintf("corr = %f, pval = %f", cor_test$estimate, cor_test$p.value), "figs/main_figures/summaries_for_paper.txt", append=TRUE)

            

## MCMC Diagnostics

all_samples %>% filter(!is.na(te), exposure_code == 1) %>% group_by(state, time, category) %>% summarize(es=coda::effectiveSize(te)) %>% pull(es) %>% quantile(., c(0.025, 0.5, 0.975))

     2.5%       50%     97.5% 
  90.5033  884.4293 1311.2619 

# all_samples %>% filter(type=="total", exposure_code == 1) %>% select(K, D, N, .chain, .iteration, te, category) %>%
#   group_by(K, D, N, category) %>% map(~ pivot_wider(.x, names_from=".chain", values_from="te"))




## Trace plots.
all_samples %>% filter(exposure_code == 1, time=="2023-01-01") %>% filter(.chain %in% c(1,2)) %>% ggplot(aes(x=.iteration, y=te)) + geom_line(aes(color=as.factor(.chain)), alpha=0.5) + facet_wrap(~category) + theme_bw()

Warning: Removed 25000 rows containing missing values or values outside the scale range
(`geom_line()`).

## Post-Hoc IMR for those who would have received abortion but did not
## 478 excess deaths
## 20,323 excess births

write_lines(sprintf("------ %i deaths over expectation divided by %i births over expectation * 1000 = %f Infant Mortality Rate for those seeking abortion who did not get one ------ ", 478, 20323, 478/20323 * 1000), 
            "figs/main_figures/summaries_for_paper.txt", append=TRUE)

### Variance Decomposition

cov_dat <- read_csv("data/dobbscovariates_2024_02_07.csv")

Rows: 4 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): state
dbl (2): median_time_post_dobbs, median_time_pre_dobbs

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ci_df <- all_samples %>% 
  filter(exposure_code == 1) %>% 
  mutate(denom = births) %>% 
  ungroup() %>% group_by(state, type, category, .draw) %>% 
  summarize(
    ypred = sum(ypred),
    outcome = sum(deaths),
    treated = sum(exp(mu_treated)), untreated = sum(exp(mu)),
    denom = sum(denom, na.rm = TRUE),
    treated_rate = treated / denom * 1000,
    untreated_rate = untreated / denom * 1000,
    outcome_rate = round(outcome / denom * 100, 2),
    outcome_diff = round(treated - untreated),
    causal_effect_ratio = treated_rate / untreated_rate
  )
ci_df_mean <- ci_df%>% group_by(type, category, state) %>% 
  summarize(causal_effect_ratio = mean(causal_effect_ratio))
ci_df_mean %>% filter(type != "total") %>% 
  group_by(type) %>% summarize(summary(lm(causal_effect_ratio ~ category))$adj.r.squared)

# A tibble: 3 × 2
  type       `summary(lm(causal_effect_ratio ~ category))$adj.r.squared`
  <chr>                                                            <dbl>
1 congenital                                                      0.0225
2 neonatal                                                       -0.0444
3 race                                                            0.199 

race_ci_mean <- ci_df_mean %>% filter(type == "race")
race_ci_mean$residual <- race_ci_mean %>% lm(causal_effect_ratio ~ category, data=.) %>% .$residual
race_ci_mean %>% group_by(state) %>% summarize(mean_residual= mean(residual)) %>% arrange(desc(mean_residual))

# A tibble: 16 × 2
   state                    mean_residual
   <chr>                            <dbl>
 1 Kentucky                      0.109   
 2 Texas                         0.0608  
 3 Alabama                       0.0283  
 4 Arkansas                      0.0212  
 5 Tennessee                     0.0157  
 6 Mississippi                   0.0151  
 7 Oklahoma                      0.0113  
 8 Louisiana                     0.00219 
 9 Ban States                    0.00211 
10 Ban States (excl. Texas)     -0.000640
11 Georgia                      -0.0160  
12 Idaho                        -0.0222  
13 Wisconsin                    -0.0306  
14 Missouri                     -0.0376  
15 West Virginia                -0.0628  
16 South Dakota                 -0.0960  

race_ci_mean %>% group_by(state) %>% summarize(mean_residual= mean(residual)) %>% arrange(desc(mean_residual)) %>% 
left_join(cov_dat, by="state") %>% select(state, mean_residual, median_time_post_dobbs) %>% ggplot(aes(x=median_time_post_dobbs, y=mean_residual)) + geom_point() + geom_label(aes(label=state)) + theme_bw(base_size=16) + ylab("Mean Residual") + xlab("Median Driving Time Post Dobbs") + geom_smooth(method="lm", se=FALSE)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 15 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_label()`).

ggsave(sprintf("figs/%ssupplement_figures/driving_time.png", fig_suffix))

Saving 7 x 5 in image
`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 15 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_label()`).

race_ci_mean %>% group_by(state) %>% summarize(mean_residual= mean(residual)) %>% arrange(desc(mean_residual)) %>% 
left_join(cov_dat, by="state") %>% select(state, mean_residual, median_time_post_dobbs) %>% ggplot(aes(x=median_time_post_dobbs, y=mean_residual)) + geom_point() + geom_label(aes(label=state)) + theme_bw(base_size=16) + ylab("Mean Residual") + xlab("Median Driving Time Post Dobbs") + geom_smooth(method="lm", se=FALSE)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 15 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_label()`).

left_join(race_ci_mean, cov_dat, by="state") %>% filter(category != "Total") %>% 
  mutate(median_time_post_dobbs = median_time_post_dobbs/60, 
         median_time_pre_dobbs = median_time_pre_dobbs/60, 
         driving_time_diff = median_time_post_dobbs - median_time_pre_dobbs) %>%
  lm(100*(causal_effect_ratio-1) ~ category + driving_time_diff, data=.) %>% summary()

Call:
lm(formula = 100 * (causal_effect_ratio - 1) ~ category + driving_time_diff, 
    data = .)

Residuals:
ALL 4 residuals are 0: no residual degrees of freedom!

Coefficients: (1 not defined because of singularities)
                           Estimate Std. Error t value Pr(>|t|)
(Intercept)                   4.966        NaN     NaN      NaN
categoryNon-Hispanic Black    9.611        NaN     NaN      NaN
categoryNon-Hispanic White    6.653        NaN     NaN      NaN
categoryOther                 5.753        NaN     NaN      NaN
driving_time_diff                NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
  (60 observations deleted due to missingness)
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 3 and 0 DF,  p-value: NA

## Compute the causal effect for southern states vs northern states
region_df <- all_samples %>% 
  filter(type == "total") %>% 
  filter(!state %in% c("Ban States", "Ban States (excl. Texas)", "Unexposed States")) %>% 
  mutate(region = ifelse(state %in% c("Wisconsin", "South Dakota", "Idaho", "Missouri"), "Non-south", "South"))

region_df <- region_df %>%
  filter(exposure_code == 1) %>% 
  ## Aggregate over south and nonsouth
  group_by(type, category, .draw, time, region) %>% 
  summarise(deaths := sum(deaths), 
                denom = sum(births), 
                ypred=sum(ypred), 
                mu = log(sum(exp(mu))),
                mu_treated = log(sum(exp(mu_treated))))

region_df %>% 
    ungroup() %>%
    ## Aggregate over time
    group_by(region, .draw) %>%
    summarize(
      ypred = sum(ypred),
      outcome = sum(deaths),
      treated = sum(exp(mu_treated)), untreated = sum(exp(mu)),
      denom = sum(denom, na.rm = TRUE),
      treated_rate = treated / denom * 1000,
      untreated_rate = untreated / denom * 1000,
      outcome_rate = round(outcome / denom * 100, 2),
      outcome_diff = round(treated - untreated),
      causal_effect_ratio = treated_rate / untreated_rate
    ) %>% pivot_wider(names_from = region, values_from = c(causal_effect_ratio), id_cols=c(.draw)) %>% 
    mutate(causal_effect_ratio_diff = `South` - `Non-south`) %>% 
    summarize(mean_rate_south = mean(`South`), mean_rate_nonsouth = mean(`Non-south`),
              mean_diff = mean(causal_effect_ratio_diff), 
              lower_diff = quantile(causal_effect_ratio_diff, 0.025), 
              upper_diff = quantile(causal_effect_ratio_diff, 0.975),
              pr_south_greater = mean(`South` > `Non-south`)) -> pr_south_greater

write_lines("------ Posterior Prob Southern vs Non-Southern States ------- ", "figs/main_figures/summaries_for_paper.txt", append=TRUE)
write_lines(capture.output(print(pr_south_greater)), "figs/main_figures/summaries_for_paper.txt", append=TRUE)

## Compute Effect Differences for NH White vs Non-White
race_df <- all_samples %>% filter(!state %in% c("Ban States", "Ban States (excl. Texas)", "Unexposed States")) %>% 
      filter(type == "race") %>% 
      mutate(black = ifelse(category == "Non-Hispanic Black", "black", "nonblack")) %>%
      filter(exposure_code == 1) %>%
      ## Aggregate over all banned states
      group_by(type, black, .draw, time) %>% 
      summarise(deaths := sum(deaths), 
                denom = sum(births), 
                ypred=sum(ypred), 
                mu = log(sum(exp(mu))),
                mu_treated = log(sum(exp(mu_treated))),
                years=mean(interval(start_date, end_date) / years(1)))

race_df %>% 
    ungroup() %>%
    ## Aggregate over time
    group_by(black, .draw) %>%
    summarize(
      ypred = sum(ypred),
      outcome = sum(deaths), years = mean(years),
      treated = sum(exp(mu_treated)), untreated = sum(exp(mu)),
      denom = sum(denom, na.rm = TRUE),
      treated_rate = treated / denom * 1000,
      untreated_rate = untreated / denom * 1000,
      outcome_rate = round(outcome / denom * 100, 2),
      outcome_diff = round(treated - untreated),
      causal_effect_ratio = treated_rate / untreated_rate
    ) %>% pivot_wider(names_from = black, values_from = c(causal_effect_ratio), id_cols=c(.draw)) %>% 
    mutate(causal_effect_ratio_diff = nonblack - black) %>% ungroup() %>%
    summarize(mean_rate_black = mean(black), mean_rate_nonblack = mean(nonblack),
              mean_diff = mean(causal_effect_ratio_diff), 
              lower_diff = quantile(causal_effect_ratio_diff, 0.025), 
              upper_diff = quantile(causal_effect_ratio_diff, 0.975),
              pr_black_greater = mean(black > nonblack)) -> pr_black_greater

write_lines("------ Posterior Prob NH Black vs Non-Black ------- ", "figs/main_figures/summaries_for_paper.txt", append=TRUE)
write_lines(capture.output(print(pr_black_greater)), "figs/main_figures/summaries_for_paper.txt", append=TRUE)

## Compute Effect Differences for Congenital vs noncongenital
con_df <- all_samples %>% filter(!state %in% c("Ban States", "Ban States (excl. Texas)", "Unexposed States")) %>% 
      filter(type == "congenital") %>% 
      filter(exposure_code == 1) %>%
      ## Aggregate over all banned states
      group_by(type, category, .draw, time) %>% 
      summarise(deaths := sum(deaths), 
                denom = sum(births), 
                ypred=sum(ypred), 
                mu = log(sum(exp(mu))),
                mu_treated = log(sum(exp(mu_treated))))

con_df %>% 
    ungroup() %>%
    ## Aggregate over time
    group_by(category, .draw) %>%
    summarize(
      ypred = sum(ypred),
      outcome = sum(deaths),
      treated = sum(exp(mu_treated)), untreated = sum(exp(mu)),
      denom = sum(denom, na.rm = TRUE),
      treated_rate = treated / denom * 1000,
      untreated_rate = untreated / denom * 1000,
      outcome_rate = round(outcome / denom * 100, 2),
      outcome_diff = round(treated - untreated),
      causal_effect_ratio = treated_rate / untreated_rate
    ) %>% pivot_wider(names_from = category, values_from = c(causal_effect_ratio), id_cols=c(.draw)) %>% 
    mutate(causal_effect_ratio_diff = Congenital - `Non-congenital`) %>% ungroup() %>%
    summarize(mean_rate_con = mean(Congenital), mean_rate_noncon = mean(`Non-congenital`),
              mean_diff = mean(causal_effect_ratio_diff), 
              lower_diff = quantile(causal_effect_ratio_diff, 0.025), 
              upper_diff = quantile(causal_effect_ratio_diff, 0.975),
              pr_con_greater = mean(Congenital > `Non-congenital`)) -> pr_con_greater

write_lines("------ Posterior Prob Congenital vs Non-Congenital ------- ", "figs/main_figures/summaries_for_paper.txt", append=TRUE)
write_lines(capture.output(print(pr_con_greater)), "figs/main_figures/summaries_for_paper.txt", append=TRUE)

## Infant mortality rate in 2019 for Ban STates
infant_mortality_rate_2019 <- all_samples %>% filter(.draw == 1, ban == 1, time %in% c("2019-07-01", "2019-01-01"), type == "total") %>% 
filter(!(state %in% c("Ban States", "Ban States (excl. Texas)"))) %>% 
select(state, mu, births) %>% mutate(expected_deaths = exp(mu)) %>% summarize(sum(expected_deaths)/sum(births)*1000) 

write_lines(sprintf("------ 2019 Infant Mortality Rate is %f ------- ", as.numeric(infant_mortality_rate_2019)), 
            "figs/main_figures/summaries_for_paper.txt", append=TRUE)